#### Blame explanation analysis, Second Lucid sample, follow-up survey ########
# Author: Ryan Kennedy

# Load needed packages
library(arm)
library(tidyverse)
library(mediation)
library(here)
library(effects)
library(patchwork)
library(stargazer)
library(ggsci)
library(ggExtra)
library(tidybayes)

# Load Data
dta3 <- read_csv(here("Data", "mediation_study2.csv")) %>%
  mutate(treatment = factor(treatment, levels = c("j", "a_a", "a_d")))

# Get general correlations
cor(cbind(dta3$jblame, dta3$s3expectations, dta3$s3whoeval, dta3$s3evaluation,
          dta3$female, dta3$dem, dta3$gop, dta3$income, dta3$white,
          dta3$ed), use = "pairwise.complete.obs")
cor(dta3$s3whoeval, dta3$s3evaluation, use = "pairwise.complete.obs")

# Run basic model to check what the relationship looks like
m1 <- lm(jblame ~ a_a + a_d + exp_trust + female + dem + gop + income + white + ed, data = dta3)
summary(m1)

# Create table of results
stargazer(m1,
          dep.var.labels = c("Blame on judge"),
          covariate.labels = c("Agrees with algorithm",
                               "Disagrees with algorithm",
                               "Trust in experts",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "main_mod_model.doc"))


# Check for heterogeneous treatment effects (CATEs) based on expert trust
m1h <- lm(jblame ~ a_a + a_d + exp_trust + I(a_a * exp_trust) + I(a_d * exp_trust) + female + dem + gop + income + white + ed, data = dta3)
summary(m1h)

# Draw 1000 simulated coefficients from the distribution in the model
m1hsim <- sim(m1h, n.sims = 1000)
# Create table of results
stargazer(m1h,
          dep.var.labels = c("Blame on judge"),
          covariate.labels = c("Agrees with algorithm",
                               "Disagrees with algorithm",
                               "Trust in experts",
                               "Agrees with algorithm * trust in experts",
                               "Disagrees with algorith * trust in experts",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "TIE_mod_model.doc"))

# Create plot of CATEs for expert trust
m1h_out <- data.frame(coef(m1hsim)) %>%
  dplyr::select(a_a, a_d, exp_trust, I.a_a...exp_trust., I.a_d...exp_trust.) %>%
  rename(Ia_a = I.a_a...exp_trust., Ia_d = I.a_d...exp_trust.) %>%
  mutate(q0a_d = a_d,
         q1a_d = a_d + (0.57 * Ia_d),
         q2a_d = a_d + (0.67 * Ia_d),
         q3a_d = a_d + (0.79 * Ia_d),
         q4a_d = a_d + (1 * Ia_d)) %>%
  dplyr::select(q0a_d, q1a_d, q2a_d, q3a_d, q4a_d) %>%
  gather(quartile, tot_effect, q0a_d:q4a_d) %>%
  mutate(quantile_name = case_when(quartile == "q0a_d" ~ "Minimum",
                                   quartile == "q1a_d" ~ "1st Quartile",
                                   quartile == "q2a_d" ~ "Median",
                                   quartile == "q3a_d" ~ "3rd Quartile",
                                   quartile == "q4a_d" ~ "Maximum"),
         quantile_name = factor(quantile_name, levels = c("Minimum",
                                                          "1st Quartile",
                                                          "Median",
                                                          "3rd Quartile",
                                                          "Maximum"))) %>%
  ggplot() +
  geom_boxplot(aes(x = quantile_name, y = tot_effect)) +
  geom_hline(yintercept = 0) +
  labs(x = "Trust in experts (Quartiles)", 
       y = "Effect of disagreeing with algorithm") +
  theme_bw()
m1h_out

# Test for moderation ATMEs using parallel regression
m1m1 <- lm(jblame ~ exp_trust + female + dem + gop + income + white + ed, data = dta3[dta3$a_a == 1,])
summary(m1m1)
m1m2 <- lm(jblame ~ exp_trust + female + dem + gop + income + white + ed, data = dta3[dta3$a_d == 1,])
summary(m1m2)
m1m3 <- lm(jblame ~ exp_trust + female + dem + gop + income + white + ed, data = dta3[dta3$a_a == 0 & dta3$a_d == 0,])
summary(m1m3)

# Calculate ATMEs from parallel regressions
a_a_coef <- coef(m1m1)[2] - coef(m1m3)[2]
a_a_se <- sqrt(summary(m1m1)$coefficients[2,2]^2 + summary(m1m3)$coefficients[2,2]^2)

a_d_coef <- coef(m1m2)[2] - coef(m1m3)[2]
a_d_se <- sqrt(summary(m1m2)$coefficients[2,2]^2 + summary(m1m3)$coefficients[2,2]^2)

# Generate tibble of results
mod_tib <- tibble(Treatment = c("Agrees with algorithm", "Disagrees with algorithm"),
                  Effect = c(a_a_coef, a_d_coef), SE = c(a_a_se, a_d_se)) %>%
  mutate(EHI = Effect + 2*SE,
         ELO = Effect - 2*SE)

# Save results for reporting in Table
write_csv(mod_tib, here("Tables", "Study 2 moderation.csv"))

# Plot of ATMEs (Figure 4D)
ATME <- ggplot(mod_tib) +
  geom_linerange(aes(x = Treatment, ymin = ELO, ymax = EHI, color = Treatment), size = 2) +
  geom_point(aes(x = Treatment, y = Effect), size = 3) +
  scale_color_npg() + theme_minimal() + theme(legend.position = "none", axis.text = element_text(size = 12)) +
  labs(x = "Treatment", y = "ATME", title = "(D) Average Treatment Moderation Effect (ATMEs)") +
  geom_hline(yintercept = 0)

# Plot regression results for control (Figure 4A)
p_c <- ggplot(dta3[dta3$a_a == 0 & dta3$a_d == 0,], aes(x = exp_trust, y = jblame)) +
  geom_point(alpha = 0.2, size = 1) +
  geom_smooth(method = "lm") +
  labs(x = "Trust in experts", y = "Blame on judge", title = "(A) Control") +
  theme_minimal()
p_c <- ggMarginal(p_c, type = "histogram")

# Plot regression results for agreeing with algorithm (Figure 4B)
p_a_a <- ggplot(dta3[dta3$a_a == 1,], aes(x = exp_trust, y = jblame)) +
  geom_point(alpha = 0.2, size = 1) +
  geom_smooth(method = "lm")+
  labs(x = "Trust in experts", y = "Blame on judge", title = "(B) Agrees") +
  theme_minimal()
p_a_a <- ggMarginal(p_a_a, type = "histogram")

# Plot regression results for disagreeing with algorithm (Figure 4C)
p_a_d <- ggplot(dta3[dta3$a_d == 1,], aes(x = exp_trust, y = jblame)) +
  geom_point(alpha = 0.2, size = 1) +
  geom_smooth(method = "lm") +
  labs(x = "Trust in experts", y = "Blame on judge", title = "(C) Disagrees") +
  theme_minimal()
p_a_d <- ggMarginal(p_a_d, type = "histogram")

# Combine plots to create Figure 4
(patchwork::wrap_elements(p_c) + patchwork::wrap_elements(p_a_a) + 
    patchwork::wrap_elements(p_a_d)) / 
  ATME

# Create table of results for SI
stargazer(m1m3,m1m1,m1m2,
          dep.var.labels = c("Blame on judge"),
          column.labels = c("Control", "Agrees", "Disagrees"),
          covariate.labels = c("Trust in experts",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "TIE_form_mod_model.doc"))

############### Mediation analysis for agreement with algorithm
# Mediation analysis for expectations in agreeing with algorithm condition
tot.fit <- lm(jblame ~ a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(tot.fit)
med.fit <- lm(s3expectations ~ a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(med.fit)
med.fit_sim <- sim(med.fit, n.sims = 1000)
med.fit_sim_df <- data.frame(coef(med.fit_sim)) %>%
  dplyr::select(a_a) %>%
  rename(a_a_med = a_a)
out.fit <- lm(jblame ~ s3expectations + a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(out.fit)
out.fit_sim <- sim(out.fit, n.sims = 1000)

# Plot of the out.fit model in bars
out.fit_sim_df <- data.frame(coef(out.fit_sim)) %>%
  dplyr::select(s3expectations, a_a) %>%
  rename(a_a_out = a_a) %>%
  bind_cols(med.fit_sim_df) %>%
  gather(treatment, sims, s3expectations:a_a_med) %>%
  group_by(treatment) %>%
  summarize(meffect = mean(sims),
            sdeffect = sd(sims)) %>%
  ungroup() %>%
  mutate(ehi = meffect + 2*sdeffect,
         elo = meffect - 2*sdeffect) %>%
  mutate(treatlabs = case_when(treatment == "s3expectations" ~ "Effect of expectations on judge blame",
                               treatment == "a_a_out" ~ "Effect of agreeing with algorithm on judge blame",
                               treatment == "a_a_med" ~ "Effect of agreeing with algorithm on expectations")) %>%
  mutate(treatlabs = factor(treatlabs, levels = c("Effect of expectations on judge blame",
                                                  "Effect of agreeing with algorithm on judge blame",
                                                  "Effect of agreeing with algorithm on expectations"))) %>%
  ggplot() +
  geom_errorbar(aes(x = treatlabs, ymin = elo, ymax = ehi)) +
  geom_point(aes(x = treatlabs, y = meffect)) + 
  geom_hline(yintercept = 0) + coord_flip() +
  labs(x = "Variable", y = "Regression coefficient") + theme_bw() 
out.fit_sim_df

# Plot effect of agreement with algorithm on expectations of accuracy
a_a_exp_plot_m <-  data.frame(coef(out.fit_sim)) %>%
  dplyr::select(s3expectations, a_a) %>%
  rename(a_a_out = a_a) %>%
  bind_cols(med.fit_sim_df) %>%
  ggplot(aes(x = a_a_med)) +
  stat_dotsinterval() + geom_vline(xintercept = 0) +
  theme_minimal() +
  labs(x = "Estimate", y = "Density", title = "(A) Effect of agreement on accuracy expectations")
a_a_exp_plot_m

# Plot effect of accuracy expectations on blame
a_a_exp_plot_o <-  data.frame(coef(out.fit_sim)) %>%
  dplyr::select(s3expectations, a_a) %>%
  rename(a_a_out = a_a) %>%
  bind_cols(med.fit_sim_df) %>%
  ggplot(aes(x = a_a_out)) +
  stat_dotsinterval() + geom_vline(xintercept = 0) +
  theme_minimal() +
  labs(x = "Estimate", y = "Density", title = "(B) Effect of accuracy expectations on blame")
a_a_exp_plot_o

# Create table of results
stargazer(med.fit, tot.fit, out.fit,
          dep.var.labels = c("Expectations","Blame on judge","Blame on judge"),
          covariate.labels = c("Expectations",
                               "Agrees with algorithm",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "aa_exp_med.doc"))

# Test mediation
med.out <- mediate(med.fit, out.fit, treat = "a_a", mediator = "s3expectations",
                   boot = TRUE, sims = 1000)
summary(med.out)

# Mediation analysis for evaluation in agreement with algorithm
tot.fit2 <- lm(jblame ~ a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(tot.fit2)
med.fit2 <- lm(s3whoeval ~ a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(med.fit2)
out.fit2 <- lm(jblame ~ s3whoeval + a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(out.fit2)

# Create table of regression results
stargazer(med.fit2, tot.fit2, out.fit2,
          dep.var.labels = c("Who evaluated","Blame on judge","Blame on judge"),
          covariate.labels = c("Who evaluated",
                               "Agrees with algorithm",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "aa_weval_med.doc"))

# Mediation analysis
med.out2 <- mediate(med.fit2, out.fit2, treat = "a_a", mediator = "s3whoeval",
                    boot = TRUE, sims = 1000)
summary(med.out2)

# Sensitivity analysis
sens.out2 <- medsens(med.out2, rho.by = 0.1, effect.type = "indirect", sims = 1000)
summary(sens.out2)
par(mfrow = c(1,2))
plot(sens.out2, sens.par = "rho", main = "Who evaluated", ylim = c(-0.2, 0.2))
plot(sens.out2, sens.par = "R2", r.type = "total", sign.prod = "positive")
par(mfrow = c(1,1))

# Mediation analysis for evaluation in agreement with algorithm
tot.fit3 <- lm(jblame ~ a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(tot.fit3)
med.fit3 <- lm(s3evaluation ~ a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(med.fit3)
out.fit3 <- lm(jblame ~ s3evaluation + a_a + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_a"),])
summary(out.fit3)

# Create table of results
stargazer(med.fit3, tot.fit3, out.fit3,
          dep.var.labels = c("Should evaluate","Blame on judge","Blame on judge"),
          covariate.labels = c("Should evaluate",
                               "Agrees with algorithm",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "aa_seval_med.doc"))
med.out3 <- mediate(med.fit3, out.fit3, treat = "a_a", mediator = "s3evaluation",
                    boot = TRUE, sims = 1000)
summary(med.out3)

# Plot effect of agreement with algorithm on desired locus of judgment
p_mfit3 <- data.frame(coef(sim(med.fit3, n.sims = 1000))) %>%
  dplyr::select(a_a) %>%
  ggplot(aes(x = a_a)) +
  stat_dotsinterval() + geom_vline(xintercept = 0) +
  theme_minimal() +
  labs(x = "Estimate", y = "Density", title = "(A) Effect of agreement on desired locus of judgment")

# Plot effect of desired locus of judgment on blame
p_ofit3 <- data.frame(coef(sim(out.fit3, n.sims = 1000))) %>%
  dplyr::select(s3evaluation) %>%
  ggplot(aes(x = s3evaluation)) +
  stat_dotsinterval() + geom_vline(xintercept = 0) +
  theme_minimal() +
  labs(x = "Estimate", y = "Density", title = "(B) Effect of desired locus of judgment on blame")

# Combine plots for presentation
p_mfit3 + p_ofit3

# Sensitivity analysis
sens.out3 <- medsens(med.out3, rho.by = 0.1, effect.type = "indirect", sims = 1000)
summary(sens.out3)
par(mfrow = c(1,2))
plot(sens.out3, sens.par = "rho", main = "Should evaluate", ylim = c(-0.2, 0.2))
plot(sens.out3, sens.par = "R2", r.type = "total", sign.prod = "positive")
par(mfrow = c(1,1))

################## Mediation analysis for disagreement with algorithm
tot.fit <- lm(jblame ~ a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(tot.fit)
med.fit <- lm(s3expectations ~ a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(med.fit)
med.fit_sim <- sim(med.fit, n.sims = 1000)
med.fit_sim_df <- data.frame(coef(med.fit_sim)) %>%
  dplyr::select(a_d) %>%
  rename(a_d_med = a_d)
out.fit <- lm(jblame ~ s3expectations + a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(out.fit)
out.fit_sim <- sim(out.fit, n.sims = 1000)

# Plot range results of analysis
out.fit_sim_df <- data.frame(coef(out.fit_sim)) %>%
  dplyr::select(s3expectations, a_d) %>%
  rename(a_d_out = a_d) %>%
  bind_cols(med.fit_sim_df) %>%
  gather(treatment, sims, s3expectations:a_d_med) %>%
  group_by(treatment) %>%
  summarize(meffect = mean(sims),
            sdeffect = sd(sims)) %>%
  ungroup() %>%
  mutate(ehi = meffect + 2*sdeffect,
         elo = meffect - 2*sdeffect) %>%
  mutate(treatlabs = case_when(treatment == "s3expectations" ~ "Effect of expectations on judge blame",
                               treatment == "a_d_out" ~ "Effect of disagreeing with algorithm on judge blame",
                               treatment == "a_d_med" ~ "Effect of disagreeing with algorithm on expectations")) %>%
  mutate(treatlabs = factor(treatlabs, levels = c("Effect of expectations on judge blame",
                                                  "Effect of disagreeing with algorithm on judge blame",
                                                  "Effect of disagreeing with algorithm on expectations"))) %>%
  ggplot() +
  geom_errorbar(aes(x = treatlabs, ymin = elo, ymax = ehi)) +
  geom_point(aes(x = treatlabs, y = meffect)) + 
  geom_hline(yintercept = 0) + coord_flip() +
  labs(x = "Variable", y = "Regression coefficient") + theme_bw() 
out.fit_sim_df

# Plot effect of disagreement with algorithm on accuracy expectations
a_d_exp_plot_m <-  data.frame(coef(out.fit_sim)) %>%
  dplyr::select(s3expectations, a_d) %>%
  rename(a_d_out = a_d) %>%
  bind_cols(med.fit_sim_df) %>%
  ggplot(aes(x = a_d_med)) +
  stat_dotsinterval() + geom_vline(xintercept = 0) +
  theme_minimal() +
  labs(x = "Estimate", y = "Density", title = "(C) Effect of disagreement on accuracy expectations")
a_d_exp_plot_m

# Plot effect of accuracy expectations on blame
a_d_exp_plot_o <-  data.frame(coef(out.fit_sim)) %>%
  dplyr::select(s3expectations, a_d) %>%
  rename(a_d_out = a_d) %>%
  bind_cols(med.fit_sim_df) %>%
  ggplot(aes(x = a_d_out)) +
  stat_dotsinterval() + geom_vline(xintercept = 0) +
  theme_minimal() +
  labs(x = "Estimate", y = "Density", title = "(D) Effect of accuracy expectations on blame")
a_d_exp_plot_o

# Combine plots for presentation
(a_a_exp_plot_m + a_a_exp_plot_o) /
  (a_d_exp_plot_m + a_d_exp_plot_o)

# Create table of results
stargazer(med.fit, tot.fit, out.fit,
          dep.var.labels = c("Expectations","Blame on judge","Blame on judge"),
          covariate.labels = c("Expectations",
                               "Disagrees with algorithm",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "ad_exp_med.doc"))

# Mediation analysis
med.out <- mediate(med.fit, out.fit, treat = "a_d", mediator = "s3expectations",
                   boot = TRUE, sims = 1000)
summary(med.out)

# Mediation analysis for preferred locus of evaluation in disagreement with algorithm
tot.fit2 <- lm(jblame ~ a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(tot.fit2)
med.fit2 <- lm(s3whoeval ~ a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(med.fit2)
out.fit2 <- lm(jblame ~ s3whoeval + a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(out.fit2)

# Create table of results
stargazer(med.fit2, tot.fit2, out.fit2,
          dep.var.labels = c("Who evaluated","Blame on judge","Blame on judge"),
          covariate.labels = c("Who evaluated",
                               "Disagrees with algorithm",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "ad_weval_med.doc"))

# Mediation analysis
med.out2 <- mediate(med.fit2, out.fit2, treat = "a_d", mediator = "s3whoeval",
                    boot = TRUE, sims = 1000)
summary(med.out2)


# Mediation analysis for preferred locus of evaluation in disagreement with algorithm
tot.fit3 <- lm(jblame ~ a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(tot.fit3)
med.fit3 <- lm(s3evaluation ~ a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(med.fit3)
out.fit3 <- lm(jblame ~ s3evaluation + a_d + female + dem + gop + income + white + ed, data = dta3[dta3$treatment %in% c("j","a_d"),])
summary(out.fit3)

# Create table of results
stargazer(med.fit3, tot.fit3, out.fit3,
          dep.var.labels = c("Should evaluate","Blame on judge","Blame on judge"),
          covariate.labels = c("Should evaluate",
                               "Disagrees with algorithm",
                               "Female",
                               "Democrat",
                               "Republican",
                               "Income",
                               "White",
                               "Education"),
          type = "html",
          out = here("Tables", "ad_seval_med.doc"))

# Mediation analysis
med.out3 <- mediate(med.fit3, out.fit3, treat = "a_d", mediator = "s3evaluation",
                    boot = TRUE, sims = 1000)
summary(med.out3)

# Sensitivity analysis
sens.out3 <- medsens(med.out3, rho.by = 0.1, effect.type = "indirect", sims = 1000)
summary(sens.out3)
par(mfrow = c(1,2))
plot(sens.out3, sens.par = "rho", main = "Should evaluate", ylim = c(-0.2, 0.2))
plot(sens.out3, sens.par = "R2", r.type = "total", sign.prod = "positive")
par(mfrow = c(1,1))

